
####-----
## Packages
####-----
rm(list=ls())
library(stargazer)
library(dplyr)
library(xtable)
library(ggplot2)
library(ggrepel)
library(cowplot) # for side by side plots in ggplot 
library(lubridate)
library(arm)


## DOWNLOAD AND SAVE DATA -- rep_data_menial.csv
## SET Working directory 
## BRING IN THE CLEAN DATA 
men_data_w_dist<-read.csv("rep_data_menial.csv")

###-----
## REPLICATION OF TABLE 3
###------

men_group <- men_data_w_dist %>%
  group_by(assembly, regionx, Total_pop.mean, typeofdistrict, chang_04_08, electricity.mean, b_type) %>% 
  summarize(num  = n())

men_group_ndc <- filter(men_group, b_type=="Pro-NDC")

r1 <-  lm(num ~ chang_04_08 + as.factor(regionx), data=men_group_ndc)
r2 <-  lm(num ~ chang_04_08 +  as.factor(typeofdistrict) + as.factor(regionx) + log(Total_pop.mean) , data=men_group_ndc) 
r3 <-  lm(num ~ chang_04_08 +  as.factor(typeofdistrict) + log(Total_pop.mean) + electricity.mean + as.factor(regionx), data=men_group_ndc)

men_group <- men_data_w_dist %>%
  group_by(assembly, regionx, Total_pop.mean, typeofdistrict, chang_04_08, electricity.mean) %>% 
  summarize(num  = n())

r4 <- lm(num ~ chang_04_08 + as.factor(regionx), data=men_group)
r5 <- lm(num ~ chang_04_08 + as.factor(typeofdistrict) + as.factor(regionx)+ log(Total_pop.mean) , data=men_group)
r6 <- lm(num ~ chang_04_08 + as.factor(typeofdistrict)   + log(Total_pop.mean) + electricity.mean +  as.factor(regionx), data=men_group)

stargazer(r4, r5, r6, r1, r2, r3,  title = "OLS regression results", dep.var.labels=c("N. of menial hires"), omit = "regionx",
          covariate.labels=c("Change in NDC vote share","Metropolitan", "Municipal", "log(Population)", "Electricity share", "Constant"), omit.stat="f", add.lines=list(c("Region fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))



####--------
### APPENDIX P ###
###---------

### Check for Pro-NPP bureaucrats 

men_group <- men_data_w_dist %>%
  group_by(assembly, regionx, Total_pop.mean, typeofdistrict, chang_04_08, electricity.mean, b_type) %>% 
  summarize(num  = n())

men_group_npp <- filter(men_group, b_type=="Pro-NPP")

r1 <-  lm(num ~ chang_04_08 + as.factor(regionx), data=men_group_npp)
r2 <-  lm(num ~ chang_04_08 +  as.factor(typeofdistrict) + as.factor(regionx) + log(Total_pop.mean) , data=men_group_npp) 
r3 <-  lm(num ~ chang_04_08 +  as.factor(typeofdistrict) + log(Total_pop.mean) + electricity.mean + as.factor(regionx), data=men_group_npp)

stargazer(r1, r2, r3,  title = "OLS regression results", dep.var.labels=c("N. of menial hires"), omit = "regionx",
          covariate.labels=c("Change in NDC vote share","Metropolitan", "Municipal", "log(Population)", "Electricity share", "Constant"), omit.stat="f", add.lines=list(c("Region fixed effects", "Yes", "Yes", "Yes")))



####--------
### APPENDIX Q ###
###---------

##-----------
## Marginal effects from reg 6
##-----------

par(mfrow=c(1,2))

holder <- matrix(NA, nrow=7, ncol = 3)

num <- c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3)
newdata2 <- data.frame(
  chang_04_08 = -0.3,
  regionx = as.factor(men_group$regionx),
  typeofdistrict= "District",  #as.factor(men_group$typeofdistrict),
  Total_pop.mean = log(men_group$Total_pop.mean),
  electricity.mean = men_group$electricity.mean) 
newdata2 <- filter(newdata2, regionx!=6)
for(i in 1:7){
newdata2$chang_04_08 <- num[i]
## Plot 
preds<- predict(r6, newdata2, type="response", se.fit=TRUE)
holder[i,1] <- mean(preds$fit, na.rm=T) # predicted
holder[i,2] <- mean(preds$fit, na.rm=T) - (1.96*mean(preds$se.fit, na.rm=T)) # lower bounds
holder[i,3] <- mean(preds$fit, na.rm=T) + (1.96*mean(preds$se.fit, na.rm=T)) # upper bounds
}


plot(c( -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,1], type="l", ylab="Number of menial hires", xlab="Change in NDC vote share (2004-2008)", bty="n", ylim=c(0, 120), main="District assemblies")
lines(c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,2], lty=2)
lines(c( -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,3], lty=2)
rug(men_group$chang_04_08[men_group$typeofdistrict=="District"], ticksize = 0.08, side = 1, lwd = 1, col = "grey50")

## Municipal 

holder <- matrix(NA, nrow=7, ncol = 3)

num <- c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3)
newdata2 <- data.frame(
  chang_04_08 = -0.3,
  regionx = as.factor(men_group$regionx),
  typeofdistrict= "Municipal",  #as.factor(men_group$typeofdistrict),
  Total_pop.mean = log(men_group$Total_pop.mean),
  electricity.mean = men_group$electricity.mean) 
newdata2 <- filter(newdata2, regionx!=6)
for(i in 1:7){
  newdata2$chang_04_08 <- num[i]
  ## Plot 
  preds<- predict(r6, newdata2, type="response", se.fit=TRUE)
  holder[i,1] <- mean(preds$fit, na.rm=T) # predicted
  holder[i,2] <- mean(preds$fit, na.rm=T) - (1.96*mean(preds$se.fit, na.rm=T)) # lower bounds
  holder[i,3] <- mean(preds$fit, na.rm=T) + (1.96*mean(preds$se.fit, na.rm=T)) # upper bounds
}


plot(c( -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,1], type="l", ylab="Number of menial hires", xlab="Change in NDC vote share (2004-2008)", bty="n", ylim=c(20, 120),  main="Municipal assemblies")
lines(c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,2], lty=2)
lines(c( -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3), holder[,3], lty=2)
rug(men_group$chang_04_08[men_group$typeofdistrict=="Municipal"], ticksize = 0.08, side = 1, lwd = 1, col = "grey50")


